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ABSTRACT 

A  second  order  scheme  for  the  solution  of  the  transient  fundamental 


semiconductor  device  equations  is  presented  which  does  not  suffer  from 
timestep  restrictions  due  to  the  stiffness  of  the  analytical  problem.  The 
second  order  accuracy  as  well  as  the  stability  properties  are  demonstrated  on 


the  simulation  of  a  p-n- junction  diode. 
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SIGNIFICANCE  AND  EXPLANATION 


In  recent  years  there  has  been  an  increasing  amount  of  interest  in  the 
calculation  of  A-C  properties  of  silicon-based  semiconductor  devices,  mainly 
due  to  the  rapid  developments  in  fabrication  technology  and  the  resulting 
miniaturization.  This  makes  the  numerical  solution  of  the  so-called  transient 
fundamental  semiconductor  equations  an  absolute  necessity.  Most  existing 
numerical  methods  for  this  problem  are  only  first-order  accurate  in  time 
because  of  the  quite  complicated  anlaytical  structure  of  the  equations 
('stiffness').  In  this  paper  we  present  a  new  method  which  is  based  on  an 
asymptotic  analysis  of  this  structure.  In  contrast  to  existing  methods,  it  is 
genuinely  of  second  order  and  does  not  suffer  from  any  timestep  restrictions. 
It  therefore  improves  the  overall  efficiency  of  A-C  calculations  significantly 
compared  to  the  existing  methods. 


A  SECOND  ORDER  DIFFERENCE  SCHEME  FOR 
TRANSIENT  SEMICONDUCTOR  DEVICE  SIMULATION 

Christian  A.  Ringhofer 


0.  Introduction 

I m  this  paper  we  present  a  finite  (iifference  discretization  method  for  the 
basic  semiconductor  device  equations  in  the  transient  case.  These  equations  - 
describing  the  field,  carrier  -  and  current  concentration  in  a  semiconductor 
material  -  are  given  by  (see  c.f.  Von  Roosbroeck  [1950]) 

(0.1)  cAY  =  q(n-p-C) 

(0.2)  a)  V  •  3  =  q[n  +R],  b)  V  •  3  =  -q[p.+R] 

n  t  p  t 

(0.1)  a)  3  =  q[D  Vn  -  ;i  nV'F],  b)  3  =  -q[D  Vp  +  u  pV'F] 

n  '  n  n  p  ^  p  p 

The  independent  variables  in  (0.1)  -  (0.3)  have  the  following  meaning: 

V  electric  potential 

-VV  electric  field 

n  negative  (electron-)  carrier  density 
p  positive  (hole-)  carrier  density 

3^  electron  current  density 

3^  hole  current  density 

t  and  q  in  (0.1)  -  (0.3)  denote  the  dielectric  permit ivity  and  the  unit  charge 

a.ud  are  material  constants.  D  .  D  and  u  ,  ii  denote  the  diffusion  coefficients 

n’  p  n  p 

.ind  mobilities  of  electrons  and  holes  which  in  general  are  functions  of  the 
electric  field  -VF.  It  is  assumed  that  (linstein’s  relation 
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holds  where  Uy  -  the  tnermal  voltage  -  is  only  dependent  on  temperature  and 
therefor  a  constant.  C(x)  In  (0.1)  denotes  the  doping  concentration  (the 
preconcentration  of  electrons  and  holes  in  the  semiconductor).  The  regions 
where  C(x)  >  0  and  C(x)  <  0  holds  are  called  the  n  and  p  regions  respectively 
and  the  manifolds  given  by  C(x)  =  0  are  called  the  p-n  junctions.  R  in  (0.2) 
denotes  the  generation  -  recombination  rate  of  electrons  and  holes  and  will  in 
general  be  a  nonlinear  function  of  n,p  and  V’F.  (See  Sze  [1969]  or  Sellcrherr 
[1984]  for  oifferent  models  for  R).  (0.1)  -  (0.3)  will  be  subject  to  the 

initial  conditions 

(0.5)  l'(x,0)  =  ’Pj(x),  n(x,0)  =  ny(x),  p(x,0)  -  p^  (x) 

where 

(0.6)  eAfj  =  nj  ■  Pj  "  C 

has  to  hold  because  of  compatibility  with  (0.1).  This  paper  is  concerned  with 
the  twodimensional  model  only.  So  (0.1)  -  (0.3)  will  hold  for  x  e  Q  C  R2.  At 
each  point  of  the  boundary  three  boundary  conditions  are  imposed  which  we 
write  in  the  form 

(0.7)  H(f ,n,p,7f,Jn,0p)  ■=  b(t)  x  e  30. 

The  shape  of  0  and  the  boundary  conditions  B  will  depend  on  the  particular 
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device  under  consideration.  A  typical  case  is  the  following: 

(0.8)  3ft  =  3ft i  J  3ft2 

(0.9)  a)  f  =  V(t),  b)  n-p~n^,  c)  n  -  p  -  C  =  0  for  x  e  3fl] 

(0.10)  a)  VTf  •  n  =  0,  b)  0  •  n  =  0,  c)  3  •  n  =  0  for  x  e  3ft? 

n  p  * 

where  3ftj  represents  Ohmic  contacts  and  3ft2  an  insulator,  n  denotes  tlie  normal 
vector  on  the  boundary  and  n^  is  a  given  constant.  V(t)  in  (0.10)  represents 
the  voltage  applied  to  the  device  which  is  varied  in  time.  (i.e.  the  input  is 
'switched' . ) 

The  scope  of  this  paper  is  to  develop  a  finite  difference  method  for  (0.1) 

-  (0.3)  which  is  second  order  accurate  in  time.  As  shown  by  Ringhofer  [1964]  b) 
the  problem  is  extremely  stiff  in  time  and  at  the  same  time  exhibits  an  internal 
layer  behavior  in  the  spaeial  directions.  This  causes  simple  one-step  methods 
based  on  central  differencing  (i.f.  Crank-Nicholson)  to  suffer  from  prohibitive 
timestep  restrictions.  First  we  seal?  the  problem  appropriately  in  Sec.  1. 

Then  -  by  maxing  use  of  a  singular  perturbation  analysis  for  the  steady  state  - 
and  the  transient  problem  carried  out  in  Markorwich  and  Ringhofer  [1984], 
Markorwich  [1984]  and  Ringhofer  [1984]  b)  we  develop  a  discretization  scheme 
which  does  not  suffer  from  any  timestep  restrictions  end  although  of  first  order 
formally  is  second  order  accurate  for  the  range  of  timesteps  one  wishes  to  use 
in  practice.  In  on  unsealed  form  it  is  of  the  form 


■  .  pi  npi  ^  J L  J  I  I  1 1 IJ I  IJI»  IJ I IJ I  IJ  ■  ■  P  ![■  M  W  U'» I  P  Vm  BP  I.  1 1. '  »  ' H 


•  ' 


w'fe'iTJ 
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(0.11)  rA»kt1  =  q[nk>1  -  Ok*'  -  C] 


(0.12) 


1  /  k +1  k, 

(n  -n  ) 


At, 


(0.13) 


1  /  k+1  k. 

(p  -p  ) 


At, 


—  V  •  [3Jk+1  +  3k  +  3k+1  -  3kj  -  1  [Rk  +  Rk+1] 
4q  n  n  P  P  2 


-  L_  [33k+1  +  Jk  +  3k+1  -  0k]  -  1  [Rk  +  Rk+1 ] 
p  p  n  n  2 


4q 


(0.14)  J?'  =  qLD^'Vn^  -  u?'n?'V'?*']  ,  9.  -  k,  k  +  1 

n  1  ii  ri  1  ' 


n  n 


(0.15)  3j  =  -  q[oV  +  i/pS/]  ,  i  =  k,  k  +  1 


whore  the  superscript  l  denotes  a  variable  at  the  discrete  timelevel  t  and 

Z  9  9, 

At.  =  t,  .  -  t.,  holds.  R  ,  D  vi  e.t.c.  denote  R,  D  ,  p  e.t.c.  with  their 
arguments  taken  at  t  =  t^.  The  spacial  differential  operators  in  (0.11  -  (0.15) 
are  discretized  by  standard  exponentially  fitted  finite  differences  due  to 
Scharfelter  and  Gummel  [1969],  The  spacial  discretization  is  given  in  Sec.  3. 

In  See.  4.  the  time  discretization  is  derived  and  an  argument  is  given  that 
(0.11)  -  (0.15)  is  second  order  accurate  as  long  as  the  timesteps  are  larger 
than  the  dielectric  relaxation  time.  The  second  order  accuracy  is  verified 
nu.r.u'  i  al  ly  i  n  Sec.  6. 


1.  Scaling  and  Singular  Perturbation  Formulation 


In  this  section  we  scale  (0.1)  -  (0.3)  appropriate ly  and  reformulate  it  as  a 


singular  perturbation  problem.  We  assume 


(1.1)  ii  -  y  =  |i  -  const 
n  p 


for  tliis  and  the  next  two  sections.  It  should  be  pointed  out  that  (1.1)  is 
probably  a  highly  unrealistic  assumption  and  therefor  inadequate.  Assumption 
(1.1)  is  only  made  for  the  sake  of  notational  simplicity  and  in  Sec.  5.  the 
scaling  and  the  scheme  is  given  for  the  case  of  a  general  model  for  the 
mobilities.  We  scale  the  independent  variables  x  and  t  to  0(1).  That  means 
that  after  scaling  the  size  of  0  and  the  speed  with  which  V(t)  in  (0.9)  is 
varied  is  0(1). 


(1.2) 


£  -  Zx  t  t  -  lz(uTy )-1t 
S  ! 


Here  l  denotes  some  characteristic  device  length.  We  scale  n  and  p  such  that 
their  maximal  size  at  Ohmic  contacts  is  0(1)  in  (0.9): 


(1.3)  n  =  Cn  ,  p  =  Cp  ,  C:=  max  C(x) 

S  3 


or  7,3  and  J  wc  use  the  sealing 
n  p 


7  =  u  7  ,  :Jn  =  qCuTp£-13n  ,  3  =  qCu^t^O  . 

s  s 


(1.4) 
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After  scaling  (0.1)  -  (0.3)  assumes  the  form 


(1.3) 


=  "s  -  Ps  '  <v  (Cs(x)  =  - ) 

max  C(x) 


( 1 .  t>)  a)n  =  V  •  0  -  R  ,  b)p  =  -  V  •  3  -R 

s.  n  s’  Ks.  p  s 

t  s  t  Ks 


(1.7)  a)  3  =  Vn  -  n  VY  ,  b)  3  =  -  V  -  p  V'F  . 

n  s  s  s  p  p  rs  s 

s  '  s  1  s 


The  subscripts  s  will  be  ommitted  from  hereon.  The  initial  and  boundary 
conditions  have  to  be  scaled  accordingly.  X  in  (1.3)  is  a  small  constant  and 
acts  as  a  singular  perturbation  parameter. 


(1.8)  X  =  (euT)1/2  (q*2Cr1/2 


For  practical  devices  X  <<  1  will  hold.  (For  instance  for  alp  silicon  device 
with  a  maximal  doping  C  =  1019  cm" 3  X  =  0(10-3)  will  hold.)  X,  which  is 
proportional  to  the  minimal  Debey  length,  lias  been  used  as  a  perturbation 
parameter  in  a  singular  perturbation  analysis  of  (1.5)  -  (1.7)  in  the  steady 
state  case  by  tiarkorwich  and  Ringhofer  [1984]  and  Markorwich  [1984]  and  for  the 
transient  problem  by  Ringhofer  L1984]  b) .  Since  we  rely  heavily  on  their 
results  for  the  development  of  our  method  we  briefly  sketch  them  in  the 
following  section. 


i 
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2.  Analytical  Results 


In  this  section  we  discuss  some  analytical  properties  of  the  solution  of 
(1.5)  -  (1.7).  These  results  will  be  used  in  the  following  sections  to  motivate 
the  choice  of  the  discretization  scheme.  We  are  motivated  by  a  singular 
perturbation  analysis  carried  out  previously  for  the  problem  (1.5)  -  (1.7).  We 
therefor  focus  on  the  information  about  the  structure  of  the  solution  obtained 
by  this  asymptotic  analysis.  For  the  precise  results  as  well  as  the  proofs  of 
these  results  we  refer  the  reader  to  the  corresponding  literature.  It  turns  out 
that  the  problem  (1.5)  -  (1.7)  is  singularly  perturbed  in  both  the  space  and  the 
time  dimension  i.e.  the  solution  exhibits  a  layer  structure  in  the  special 
variables  arid  can  evolve  on  several  different  timescales  simultaneously.  Away 
from  t  =  0  and  the  layer  regions  the  solution  of  (1.5)  -  (1.7)  will  be  given  up 
to  terms  of  order  0(a)  by  the  outer  solution  w  -  (  T,n ,  p,^  ,5^)^  of 


(2.1)  0  =  n  -  p  -  C 


(2.2)  a)  n  =  V«3  -  It  ,  b)  p_  =  7*3  -  R 

t  n  rt  p 


(2.3)  a)  3  =  Vn  -  nvy  ,  b)  J  =  -Vp  -  pVf. 

n  ’  p 


So  Poisson's  equation  degenerates  for  A  ->  0  into  the  condition  of  vanishing 
space  charge.  Because  of  this  drop  in  order  w  can  not  be  expected  to  satisfy 
all  boundary  conditions  (0.7)  automat ical ly.  It  turns  out  however  that  (2.1)  is 
consistent  with  (0.9)  c)  and  one  of  the  conditions  (0.10)  a),  b)  is  satisfied 
automatically  as  soon  as  the  other  one  is  because  VC*n  =  0  holds  on  the 
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oppropriate  parts  of  the  boundary  3ft.  Therefor  ail  boundary  conditions  (0.9)  - 
(0.10)  can  be  satisfied  by  w  and  no  boundary  layers  occur  (see  Markowich 
[l9o4]).  We  are  confronted  however  with  an  internal  layer  structure  around  the 
p-n  junctions.  The  doping  profile  C(x)  will  vary  very  rapidly  across  the 
p-n- junct ions  and  in  the  extreme  case  is  even  modelled  as  discontinuous  there. 

In  either  case  n  and  p  have  to  be  discontinuous  across  the  p-n- junctions  in  order 
to  satisfy  (2.1)  and  this  causes  the  'full'  solution  w  to  develop  a  layer 
there,  because  w  is  discontinuous  (2.2)  -  (2.3)  can  not  hold  at  the 
p-n- junct  ion  and  the  reduced  probLem  has  to  be  supplemented  by  4  jump  conditions 
there.  These  conditions  have  been  derived  in  Markowich  and  Ringhofer  [1984]  and 
Markowich,  Ringhofer  et.  al.  [1983]  and  are  of  the  form: 


(2.4) 


-  'F  4*  . ,  >  n  ♦ 

ne  ,  pc  ,  Jn*  n,  J  •  n 


are  continuous  across  the  p-n-junction  where  n  denotes  the  unit  normal  vector  on 
the  p-n-junction.  Ihese  conditions  say  that  the  currents  and  the 
ct|ui  -  fc  rmi- paten  t  ials  are  preserved  across  the  junctions.  Inside  the 
p-n-junction  la;,  ers  7  n  and  p  vary  in  the  ’fast'  spacial  variable  £;  =  d/  \  where 
<i  denotes  the  perpend  i  cu  I  or  distance  of  a  point  x  from  the  p-n-junction.  They 
i r i •  tjivn  by  n,  p  satisfying 


Us,t)<- 


(2.4) 


i )  n 


b)  p  -  B(s,t)e 
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where  s  is  some  parametrization  of  the  p-n-junction  and  A  and  U  are  determined 

bv  the  solution  w  of  the  reduced  probLem  (2.1)  -  (2.4).  3  and  3  do  not 

n  p 

exhibit  any  layer  behavior.  The  spacial  structure  of  the  solution  is  depicted 
graphically  in  the  figure  below. 


p-n-junction 


n  =  A  ( s ,  t )  e  ** 
p  =  3(s,t)e 


bo  far  we  have  described  the  structure  of  the  solution  away  from  t  =  0.  If  the 
described  asymptotic  structure  is  to  hold  for  all  t  >  0,  the  initial  functions 
?^,nL  and  p(  apparently  have  to  satisfy  the  relations  (2.1)  and  (2.6).  Moreover 
dif ferentiating  (2.1)  with  respect  to  t  and  using  (2.2)  a)  b)  yields 

(3.7)  7  •  (3  43  )  =  0 

n  P 

which  gives  an  additional  condition  on  ?  ,  nj.  and  pj.  in  ilarkowich  and 
H i  nghofer  [193!ij  and  Markowich  [1934]  it  is  shown  that  these  conditions  arc 
satisfied  if  7(  ,  n^,  are  the  solutions  of  the  steady  state1  problem 


( x(s) ,y(s) )  •  d  = 


(x, 


V  =  n 

KK 

3=3 
n  n 

3=3 
P  P 


-  P 


w(x,y) 
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6.  Numerical  Results 


In  this  section  we  verify  the  stability  and  convergence  properties  of  the 
discussed  scheme  for  the  case  of  a  p-n- junction  diode.  The  length  and  depth  of 
the  device  are  taken  to  be  Ip  and  0.2b  p  respectively.  The  doping  profile  C  is 
taken  to  the  constant  10in  cm- 3  in  the  n-reyion  and  -It)1  ’  cm-1  in  the  p-reqion. 
(So  v\e  consider  an  abrupt  p-n- junction) .  After  employing  the  scaling  of  Sec.  1 
we  obtain  the  system  (l.b)  -  (1.7)  posed  on  a  rectangular  ooinain  of  the  form 


/ 

/ 

/ 

/ 

/  Insulator 

/ 

/ 

/ 

/ 

/ 


The  length  and  width  of  the  domain  are  1  and  0.2b  respectively.  X  turns  out  to 

be  b.ll)-3  in  this  case.  ihe  contacts  correspond  to  the  part  of  the 

boundary.  So  the  boundary  conditions  (u.9)  are  imposed  there.  Tbs  insolator 

corresponds  to  ri,K-*  ('J.10)  is  imposed  on  this  part  of  the  boundary.  The 

voltage  applied  to  the  diode  is  given  by  the  difference  in  f  between  the  two 

contacts  ca  I  i  hr  a  t.  •  d  by  i.lm  bu  i  1 1  in  potential  y,  (see  c.  f.  S/e  11969)).  So  the 

i 

ooutidary  condition  for  i'  in  (U.9)  a)  assumes  the  form 


(6.1)  a)  . i  (\)  ■  J-u(t)  \  d  the  .inode 

t  >  4  1 


Cathode  Contact 


Anode  Contact 
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n  -  r,  -  C 


n  V-:)  +  H  b)  |)  s  -7«3  -  R 

t  n  1  l  p 


')  =  L)S(7n  -  r i V  !' )  b)  3  -  -DS(Vp  pvy) 

n  n  p  p 

(x  ,7'?  )  =  (17  17 )~1/‘~  urU*  ,  -I  V'P  )  ,  r  =  n,p. 

o  j  1 1  p  r  5  £  5 


( 5 .  •+ )  -  (5.6)  into  the  discret  i/ation  method  (4.13)  - 
17)  has  t.o  be  replaced  by 


ir(  <’J.  n 
n  h 


( :  ln)  '7 h  3 


-b.’i  .  f  (Mp)V.  t'l 


b°  :  I)  (x.  .,7  7..)  ,  r  -  n,p 

r  r  lj’  h  ij 


(4.17)  , 
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5_.  Vjir  iabl  c  Mobilities 

\s  nun  t  Lone']  in  Sec.  1  the  assumpt.  ion  g  -  n  =  const  and  D  =  D  =  const 

n  p  n  p 

is  quite  unreal  i;t  ic.  As  [minted  out  in  Mock  [1981]  it  might  sometimes  be 
just  i  Cy.it) Le  in  steady  state  simulations  while  it  is  probably  inadequate  in  the 
transient  ease,  it  is  however  reasonable  to  assume  that  Einstein's  relation 


(8.1) 


I) 

n 


holds  where  u-j  only  depends  on  temperature  (and  therefore  is  constant  for  our 
purposes).  Various  models  for  p^  and  p  have  been  proposed  (see  c.f.  Selberherr 
[1934]  for  an  overview).  In  the  most  general  case  g^  and  g  will  depend  on  x 
and  on  the  electric  field  -  7“  i.e. 
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(4.10)  in  the  given  form  since  it  involves  a  fourth  order  difference  operator 

in  space  and  it  is  difficult  to  implement  nonconstant  mobilities  ii^  and  Up  in 

this  form.  We  therefor  rewrite  (4.7)  -  (4.10)  using  the  original  variable  set 

('I'.n.p.J  )^  and  obtain 

n  p 


(4.13) 


divhVh 


.k  +  1 


k+1 

n 


k+1 

P 


-  C 


1  ,  k+1  k.  1  r,.  .^k+l  nk  -,k+1  ..k,  ,uk  0;.k+1  n 

(4.14)  -  (n  -n  )  =  —  Ldiv, (33  +  3  +  J  +  3  )  -  2R  -  2R  J 

4  h  n  n  p  p 


/,  * .  .  1  /  k  +  1  k.  1  ,  /ank+l  nk  nk+1  nk.  ,Qk  •jDk+^i 

(4.1b)  -  (p  -p  )  =  -  —  Ldiv. (33  +3+3  -  3  )  +  2R  +  2R  J 

,',V  4  n  p  p  n  n 


0,0  0  0  0 
(4.1b)  r  Lx"7hn  '  -  (Mn ")7(  V ] 

(4.17)  3p  =  L<?'vhpp'  +  (mpr)  vhyeJ 


l  -  k,  k+1 


We  have  already  argued  heurist ical ly  that  the  method  is  second  order  accurate  in 

time  as  long  as  A2  <<  At^  holds.  This  will  be  verified  numerically  in  Sec.  6. 

for  actual  computations  (4.16)  and  (4.17)  are  inserted  into  (4.14)  and  (4.1b) 

k  k 

and  3  and  3  are  computed  after  the  fact.  Given  7  ,  n  arid  p  solving  (4.13)  - 

ri  (j  ' 

k  +  1  k+1  k  +  1 

(4.17)  amounts  to  solving  3  coupled  I’.D.L.'s  for  7  ,  n"  and  p  at  each 

time  step  which  is  the  same  arnmounl  of  labor  as  is  reguired  by  backward 
differences  at.  each  t.ime.step. 
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idea:  'Ac  use  backward  differences  for  (4.2)  a)  and  centered  differences  (the 


trapezoidal  rule)  for  (4.2)  b).  So  we  choose  a  mesh  0  =  tg  <  tj  <  ...  and 


meshsizes  6t^:-  -  t^,  k  =  0,1 

(y ,u , a , 3U)  at  t  =  tk  by  (wpj^dS 


....  and  approximate  the  gr  idf unct ions 
(4.2)  a)  -  b)  is  now  discretized  by: 


2 

(4.7)  _L  div.V.  (7k+1  -  4,k)  =  div  dk+^ 

Atk  h  h  h 

1  /  k+1  k.  1  ,  ,nk  ..k+lvi  t,k  ,,k+1 

(4.6)  -  (u  -  u  )  =  —  Ldxv.  (d  +  J  )J  -  f(  -  R 

At  2  h  u  u 

k 

(4.9)  dk  +  1  =  -(Mukf1)Vh4-k  +  1  +  <k  +  1['7h(C+X2  div^T)] 

(4.10)  d*  =  -  M(C  +  X2  div.7 u'i’a)  •  V.re 

u  h  h  h  h 

l  -  k,  k  +  1  where 

(4.11)  Rl  =  R(?h'/',n\p?'),  k'  =  <(7^*),  9.  =  k,  k  +  1 


holds.  The  local  discretization  error  L  in  (4.7)  and  (4.8)  is  given  by 

(4.12)  L  =  ( X2 At.  div  7,  f.  .  ,  At2u  ) 

k  h  n  tt  k  ttt 


Thus  uL'l  <  const  (X2At,  +  At2)  holds.  Since  we  aim  to  use  stepsizes  At,  >>  X2 

K  K  K 

L  n  -  0 ( At 2 )  holds  for  all  pract  ical  purposes.  (This  is  verified  numerically 
in  Sec.  6.)  lor  numerical  purposes  it  is  inconvenient  to  actually  solve  (4.7) 
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n.  o.  J  ,J  can  bo  recovered  bv  virtue  of 

e  p 

n  -  —  div.  V  7  +■  C)  ,  P  -  —  (u-A2  div  V,  7  -  C) 

2  h  h  ?  h  h 

14  .  J  i 

J  =1(3+3)  ,  3  =1(3-3) 

n  2  u  p  2  u 

the  bouiuidry  conditions  (3.17)  have  to  be  transformed  appropriately.  In  this 
form  it  is  clear  that  (4.2)  a)  together  with  (4.3)  constitutes  the  fast  time 
scale  confxment  while  (4.2)  b)  together  with  (4.4)  constitutes  the  slow  time 
scale  component.  Inserting  (4.3)  into  (4.2)  a)  gives  in  leading  order 

(4.6)  A2  div.  V,  7.  =  -div.  [  (Mu)  V.  '?]  +  .... 

1 1  h  t  h  h 

which  is  the  discrete  equivalent  of  (2.12).  Since  u  =  n  +  p  will  be  positive 
and  bounded  away  from  zee o  throughout  most  of  the  device  7  can  evolve  on  the 
time  scale  t  =  t/A2. 

On  the  other  hand  u  apparently  evolves  on  the  ('slow')  t-time-scale  only 
(which  also  is  the  scale  on  which  the  applied  voltage  V  in  (0.9)  is  varied.  V.c 
‘in  refer  employ  a  strategy  which  has  been  successfully  employed  by  Kreiss  and 
Oii-hois  ,197'jj  and  dinqiiofer  |  I9S4J  a).  Hint  is  to  use  lower  order,  onesided, 

! . -  •. t  lhie  schemes  for  the  fast  components  of  Lin'  system  and  to  use  symmetric 

-if.  for  the  'slow'  ioi!i|x)fi(nL  s  which  art-  of  higher  order  but  only  A-s  tabic. 

! >  in  in  general  be  expected,  that  the  strong  stability  properties  of  the 
,  is!'  m»-t  hod  can  he  ret  .lined  while  i  Is  lower  order  accuracy  does  not  destroy 

1  to-  ivriM  neeijraey.  y,e  restrict  our1. elves  to  the  simplest,  version  of  this 
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4.  Discretization  in  Time 

lie  non  der  ive  <i  discretization  scheme  for  the  system  of  ordinary 
differential  equations  (3.14)  -  (3.16).  It  turns  out  that  second  oruer  methods 
based  on  centered  differencing  (i.e.  the  trapezoidal  rule  or  the  box-scheme) 
experience  severe  stability  problems.  This  is  due  to  the  presence  of  several 
different  time  scales  in  the  analytical  problem  (1.5)  -  (1.7)  (see  Sec.  2).  So 
the  O.D.F.  system  (3.14)  -  (3.16)  is  stiff.  In  order  to  circumvent  this  problem 
(3.14)  -  (3.16)  is  transformed  into  an  equivalent  system  which  allows  the 
separation  into  a  'fast'  and  a  'slow'  component.  For  the  fast  component  we  use 
an  L-stable  Ist  order  method  (i.e.  backward  differences).  For  the  slow 
component  we  use  an  A-stable  2nc*  order  method  (i.e.  the  trapezoidal  rule).  This 
will  avoid  any  restrictions  on  the  timesteps.  On  the  other  hand  it  is  well 
known  that  in  such  a  situation  the  use  of  a  lower  order  method  for  the  fast 
component  does  not  destroy  the  over  all  (second  order)  accuracy  of  the  scheme 
precisely  because  of  the  stiffness  (see  c.f.  Kinghofer  [i?S4j  a).  So  we  obtain 
a  second  order  method  which  does  not  suffer  from  prohibitive  timestep 
restrictions.  In  order  to  separate  the  fast  and  the  slow  components  of  (3.14)  - 
(3.16)  we  differentiate  (3.14)  with  respect  to  t  and  introduce  the  new  variables 


(3.  I.S) 


(d)  3  =  [<7,  n  -  (Mn)  V.  f  J 

n  h  h 

(3.1b)  (b)  Jp  =  -  [<Vhp  +  (Mp)VhY] 

for  (x.,y^ )  e  Gj 

(3.17)  6('F,n,p,VhY,J|i,0p)  =  0  for  x±, y^)  e  G^. 

Thus  we  have  replaced  the  mixed  parabolic-elliptic  system  (1.5)  - 
ODE  system  (3.15)  -  (3.16)  coupled  to  the  alyebraic  system  (3.14). 
Lite  next  section  will  be  to  discretize  (3.14)  -  (3.16)  in  time  by 


(1.7)  by  the 
The  scope  of 
L-stable 


second  order  finite  difference  scheme. 
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The  discretization  (3.8)  -  (3.9)  can  be  motivated  in  the  following  way: 

For  a  certain  set  of  boundary  conditions  (corresponding  to  the  equilibrium 

case)  3  and  3  will  vanish  identically  and  n  and  p  will  be  of  the  form 
n  p  J  1 

(3.13)  n  =  Ae^  ,  p  =  Ae-^ 

for  some  constant  A.  In  this  case  (1.7)  a),  b)  is  solved  exactly  by  (3.8)  - 
(3.9).  In  Markowich  Ringhofer  et  al  [1983]  a  singular  perturbation  analysis  of 
the  steady  state  problem  has  been  carried  out  and  it  has  been  shown  that  <  in 
(3.11)  is  precisely  the  correct  fitting  factor  if  (1.5)  -  (1.7)  is  viewed  in  the 
framework  of  exponentially  fitted  schemes  for  singularly  perturbed  problems  (see 
c.f.  Doolan  et.al.  [1980].  (3.6)  -  (3.7)  together  with  (3.8)  -  (3.9)  together 

with  the  boundary  conditions  (0.7)  now  give  the  discretization  of  (1.5)  -  (1.7): 

(3.14)  X2  div.  7.  V  -  n  -  p  -  C 

h  h 

(3.15)  a)  —  n  --  div,  3  +  R  ,  h)  i-  p  -  -div,  3  +  R 

3t  >»  n  DL  h  R 


discretized  by  centered  differences  in  a  straightforward  manner 


(3.6)  X2  div.  7.  V  -  n  -  p  -  C 

n  h  r 

(3.7)  a)  —  n  =  div  3  -  R  ,  b)  —  p  =  -di  v.  3  -R 

St  h  n  St  h  P 

(x.,yj  c  Gj 

where  R..  =  R(n..p..)  Isolds.  ft  is  the  current  relations  (1.7)  a),  b)  which 
LJ  ij  ij 

have  to  be  discretized  with  rare.  As  mentioned  in  the  previous  section  layers 
can  be  expected  to  occur  in  ’f,  n  and  p  across  tne  p-n-junction  (where  C(x,y) 
c.h  anges  its  sign  rapidly).  Mnrkowich  and  Ringhofer  et.al.  [  1 9  S  3  ]  [joint  out  that 
a  mesh  which  eguidistr ibutes  the  local  discretization  error  will  inevitably 
contain  a  prohibitive  arninount  of  [joints  if  the  p-n-junction  does  not  lie 
parallel  to  the  gridline  directions  (which  can  not  be  expected  in  general). 

’  enarfetter  and  Gummel  [1969]  suggested  the  following  discretization  of  (1.7) 
a) ,  b) : 

'  i .  h )  J  L  r7.  n  -  (! in )  •  7  o] 

n  h  h 


r  '!  nr  *1  is  ill  •  f  i  ih  I  1)V 
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(3.3)  3  =  ^X,3y)T  ,  3X(J: - 111  ,  y.,t)  -  3X 

S  S  S  5X  2  ^  j  7  '  5 


i+1/2,  j 


3y(xt,  _i — —111)  -  3y  (t),  s  =  n,p. 


i,  j  +  1/2 


In  order  to  write  the  discrete  approximations  of  the  spaciai  differential 
operators  in  (1.5)  -  (1.7)  in  a  convenient  form  we  introduce  the  following 
difference  operators:  for  a  vectorfield  v  given  on  the  intermediate  points  of 
the  grid  G  we  define  the  discrete  divergence  operator  by 


(3.4)  d)vh  Vj.:=  2<{xt»«xl_1)  + 


2UY!S'j-1)‘1<vl,>1/2  '  ’',’erC 


v.  .  -  (vX  .  vy  .  .  ,,)  holds 

lj  1+1/2, j’  i» J+1/2 


for  a  scalar  gridf unction  u  we  define  the  discrete  gradiend  7^  by 


(3.5)  Vij:r  (d|,ui  +  l/2,j’  dhui,j+1/2)  » 


X  -  -  I  /  V 

o.  u .  .  :  6x.  (u.  +  . -u  .  . ) 

h  1+1/2, j  l  1+1, J  lj 


dhUi,j.1/2:‘  5yj  <Ul,J»1'UiJ)' 


Poisson's  e<|uation  (1.5)  and  the  continuity  equations  (1.6)  a),  b)  are  now 
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3.  Spacial  Uiscret i/at ion 


In  this  section  we  discuss  the  discret  Lz.it  ion  of  (1.5)  -  (1.7)  in  the 
spacial  variable  x  leaving  the  time  t  continuous.  Thus  (1.5)  -  (1.7)  is 
replaced  by  a  large  sparse  system  of  ordinary  differential  equations.  We 
utilize  an  exponentially  fitted  finite  difference  scheme  due  to  Scharfetter  and 
Hummel  [1 969]  which  has  been  widely  used  for  the  solution  of  the  steady  state 
problem  (2.8)  -  (2.10).  first  we  choose  a  nonuniform  but  rectangular  mesh  in 
space: 


(3.1)  C  -  {(x.,yj),  i  -•  0(1  )N ,  j  =  0(1)H}  wi' 


0  ~  XQ  <  X-]  <••••<  XN  -  1 


0  =  y0  <  Vi  <  —  <  yM  =  d 


6x  x  -  x  ,  i  =  0 ( 1 ) N - 1 ,  <ty  -  y  -y  ,  j=0(1)M-  1 
i  i+1  i  j  j+1  j 


For  obvious  reasons  it  is  convenient  to  distinguish  between  boundary  and 
interior  points.  Go  we  define 


(3.2)  G,:.  |x.,y.),  1-KDN-l,  j  =  1(1)11-1},  CQ  =  G  - 


'i ,  n  and  p  are  approximated  by  their  values  at  the  gridpoints  (x.,y^) 


w  (  x  ,  y  t)  -  w  (t)  v.  -  t' , n , p , 

J  1 J 


th.,*  current,  densities  .aid  3  are  approx  i  rated  at  the  midpoints  between  two 
gr  i  dj/> ints 
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in  leading  order 

(2.12)  w  -  -  (n+p)A7  +  X2A27  +  Jl.o.t. 
t  =  t/X2 

In  Kinghofer  [1984]  b)  a  complete  asymptotic  expansion  of  the  solution  w  on  the 
fast  timescale  t  =  t/X2  is  given  in  the  case  of  1  space  dimension.  It  turns  out 
that  only  f  and  the  drift  currents  -nV7,  -pV7  can  evolve  on  the  fast  timescale 
while  n  and  p  can  only  evolve  on  the  ’slow'  t-timescale. 
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(2.8)  X2AF=n-p-C 


(2.9)  a)  7-3  =  R  ,  b)  7*3  =  -R 

n  p 

(2.10)  a)  3  =  Vn  -  nVV  ,  b)  7*3  =  -Vp  -  pV’F 

n  p  ' 


In  Kincjhofer  [1784]  b)  it  is  shown  that  each  time  (2.7)  is  violated  the  solution 
w  wili  evolve  on  an  additional  'fast*  timescale  t / X 2  which  corresponds  to  the 
dielectric  relaxation  time  in  unsealed  form.  This  is  of  importance  for  the 
numerical  solution  of  (1.5)  -  (1.7)  because  it  means  that  the  problem  is  stiff 
in  time  and  special  stability  (i.e.  I-  -  stability)  properties  will  be  required 
of  discretization  methods.  In  Ringhofer  [1984]  b)  an  asymptotic  expansion  for 
the  solution  on  all  timescales  and  for  general  initial  functions  has  been 
carried  out.  We  only  briefly  sketch  the  results  here:  It  is  beneficial  Lo 
transform  the  system  (1.5)  -  (1.7):  Differentiating  (1.5)  with  respect  to  time 
yields 


(2.11)  X  ^  A  'F ,  =  7.(3  +  3  ) 

t  n  p 


(2.11)  states  the  conservation  of  the  total  current  consisting  of  the  drift  and 

diffusion  current  for  electrons  and  holes  (3  e3  )  and  of  the  displacement 

current  -A'V'F  .  It  is  well  known  that  the  displacement  current  can  evolve  on 

the  timescale  of  the  diclcctic  relaxation  Lime.  Whenever  7* [3  +3  ]  t  0 ( X 2 ) 

n  p 

holds  away  from  the  p-n- junction,  where  A’F  is  hounded  for  X->0,  T  will  evolve  on 
the  timescale  t / X 2 .  inserting  the  expressions  for  3  and  3  from  (1.7)  gives 


<6.])  I))  7  ( x ,  t )  =  (x)  x  at  the  cathode 

i  ^ 


We  simulate  a  switch  from  u  _  0  (equilibrium)  Lo  u  -  -0.3V  (forward  bias)  in  a 
t  imoi ntorva i  of  10~r)  see  which  corresponds  to  0.1  in  our  scaling,  for  the 
recombination  rate  the  ShockLy  Head  Hail  -  recombination  term 


(»>.?) 


n.p  -  n 


r  Um  n .  )  *-  t  (|i*ii.) 
n  r  ,)  l 


has  been  used.  (n^,  and  r  are  given  constants.)  The  nonlinear  system  of 

equations  for  (  v  ,n"  , )  arising  in  (4.13)  -  (4.17)  at  each  timestep  has 
been  solved  by  llewton's  method  using  a  damping  strategy  due  to  Deuflhard 
L  I '3 7 4 J .  Ihe  purpose  of  this  testcxample  is  twofold:  firstly  it  should  verify 
the  stability  properties  mentioned  in  Sec.  4,  i.e.  we  should  observe  that  no 
oscillations  occur  if  At^  >>  X2  is  chosen  at  any  point  in  time.  Secondly  it 
should  verify  the  second  order  accuracy  in  time.  Therefor  a  fixed  mesh  is 
chosen  in  the  spacial  direction  and  it  is  verified  numerically  how  well  the 
system  of  ordinary  differential  equations  (3.14)  -  (3.16)  is  approximated  by 
(4.13)  -  (4.17)  by  varying  the  tiinesteps  At^.  So  we  choose  a  rectangluar 
nonun i f orm  gr  id 


'•  :  f  (x.  ,y  J  ,  i  =  0(  1  )n  ,  j  =  U(  1  )m} 


where  tin-  gridlines  are  clustered  around  the  p-n- junet  ion.  The  maximal  mesh 


spat:  i  ng  h  <•'.  been  chosen  as 


(6  ) 


6x :  -  inix  1 6x .  :  6x  =  xi  +  l  '  *  >  i  =  0(1)n  -  1}  =  0.05 


6y:---  max  { 6y ^  :  6y^  =  y^+1  -  y ^ ,  j  =  0(  1  )m  -  1 }  =  0.05 


Then  a  solution  is  computed  on  a  quite  course  mesh  in  time  which  is  obtained  by 
bounding  the  yrowth  in  the  drift  -  and  diffusion  current  density:  i.e. 


(6.5) 


!<  + 1  k  k 

id  -  1!  2  <  WP0K|l2 


J*  =  3Z  +  3*  ,  i  =  k,  k  +  1 

n  p  ’  ’ 


is  required  at  each  step  with  some  suitable  parameter  w.  (Here  It  •  ll 2  denotes  the 
Discrete  L2-nurn.)  Using  this  strategy  the  coarsest  mesh 


(6.6)  T0:=  {tj  :  i  =  0(1 )K} 


is  obtained,  lhen  the  meshspacings  are  halved  successively  to  obtain  the  meshes 


l  >  1  2 


2*K}  ,  Z  = 


i  =  2j 


+  tj  )  ,  i  =  2j  +  1 


Then  the  solution  on  the  finest  mesh  is  regarded  as  'exact'  solut  ion  and  the 


relat  ive 

differences  (in  the  discrete 

1-2  norm)  to 

the  coarser- 

mesh  solutions  . 

o  omp  a  t.  ed . 

If  the  met hod  is 

of  second 

order  they 

should  decay 

asymptotically 

a  factor 

4  from  T  ^  to  J  }+y 

Table  1 

gives  the  relative  errors 

in  the  L?  norm 

TABLE  1 

l 

ESI 

N 

P 

J  =  3  +3 

n 

3P 

P 

0 

U.33E-01 

o.«yr+oo 

0. 1 81  +00 

0. 88L+00 

0.39E+00 

1 

0.10L-01 

0.17!  oo 

0.461-01 

0.1 81. -00 

0.1 21'^  00 

p 

0. 16L-U2 

0.27L-01 

0.74L-02 

0.29E-01 

0.20E-01 

3 

0.17E-03 

0.62E-02 

0.77E-03 

0.31E-02 

0.20E-02 

4 

0.39E-04 

0.  I6E-02 

0.18E-03 

0.80E-03 

0.45E-03 

8 

0.34E-08 

0.33C-03 

0.3 BE -04 

0 .2  IE-03 

0.94E-04 

As  a  test  (3.14)  -  (3.16)  has  been  solved  on  the  same  meshes  in  the  spacial  and 
the  time  direction  by  usinq  the  same  spacial  discretization  but  simple  backward 
differences  for  n^  and  pL  in  ( 3 . 1 1> )  a)  b).  The  relative  differences  to  the 
f.i  nost-gr  id-so  lut  ions  for  this  case  should  decay  by  a  factor  2  from  to  T^ 
because  of  the  first  order  accuracy.  They  are  given  in  Table  2: 


TABLE  2 


PS  I 

N 

p 

3  =  3  +3 

31  p 

3P 

u 

0.  1  li  -01 

0. 1  31 

+01 

0.821  -01 

0.701  +00 

0.14L+00 

} 

U  .6/1  -0,’ 

U .  7  8! . 

+00 

0. 301  -01 

0.171  <00 

(,.8  1!  -0  1 

4 

0. 171  -02 

0.  19!. 

•  00 

0.17L-01 

0.64E-01 

0.48L-01 

i 

0. 1-  i  -02 

0 .  !  9! 

■00 

0.831  -07 

0.  171  -01 

0 .7 3i  -0  1 

4 

0.6  !!  -01 

0.811 

-01 

0.371  -07 

0.141  -01 

0. 101  -0  1 

x> 

0.7  SI  -0  ) 

0.78! 

-01 

0.171-0? 

0.46!  -07 

0.341  -07 
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Comparing  Tables  1  and  2  it  can  also  be  seen  that  to  compute  c.f.  .7  -  3  *  3 

n  p 

with  a  relative  accuracy  of  1%  4  times  as  many  timesteps  are  required  by 
backward  differences  shown  by  the  proposed  method.  In  Figures  1  -  9  plots  of 
the  solution  at  various  time  levels.  The  switch  from  u  =  O.v  to  v  0.5v  is 
performed  in  the  interval  from  t  -  0.1  to  t  =  0.2. 


Figure  1  shows  the  initial  solution  4' ( x , 0 )  which  has  been  taken  as  the 
steady  state  solution  for  u  =  O.v  (equilibrium). 

Figures  2-8  show  the  solution  at  t  =  0.2  (the  end  of  the  switching  interval). 
In  particular  they  show 


Figure  2:  7(x,0.2) 

Figure  3;  n(x,0.2) 

Figure  4;  p ( x ,  0.2) 

Figures  5  and  6:  The  x  and  y  components  of 
Jn(x,0.2) 

Figures  7  and  8:  the  x  and  y  components  of 

0  (x ,0.2 ) 

P 

Figure  9  shows  the  total  current  flow  through  the  cathode  contact  as  a 
function  of  time: 


JT0tU) 


/  (-77  J 

cathode  ' 


3  ) 


P 


> 

n 


(is 


whore  n  is  the  normal  vector  to  the  cathode. 
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